PFG 5005 - Estudo Dirigido 1: Seções de Poincaré para a Hamiltoniana de Hénon-Heiles¶

Neste notebook, implementaremos um método de integração numérica simplético para a Hamiltoniana de Hénon-Heiles e o algoritmo de Hénon para construir seções de Poincaré.

5.1. Equações do Método de Euler Simplético¶

A Hamiltoniana de Hénon-Heiles é dada por:

$$H=\frac{1}{2}(p_{1}^{2}+p_{2}^{2}+q_{1}^{2}+q_{2}^{2})+q_{1}^{2}q_{2}-\frac{1}{3}q_{2}^{3}.$$

Para aplicar o método de Euler simplético, precisamos das equações de movimento, que são as derivadas da Hamiltoniana em relação às posições e momentos.

$$\dot{q_1} = \frac{\partial H}{\partial p_1} = p_1$$ $$\dot{q_2} = \frac{\partial H}{\partial p_2} = p_2$$ $$\dot{p_1} = -\frac{\partial H}{\partial q_1} = -(q_1 + 2q_1q_2)$$ $$\dot{p_2} = -\frac{\partial H}{\partial q_2} = -(q_2 + q_1^2 - q_2^2)$$

In [6]:
import numpy as np
import matplotlib.pyplot as plt
import os
import warnings

# converte todos os RuntimeWarning do numpy em exceções
warnings.filterwarnings("ignore", category=RuntimeWarning)
np.seterr(over='raise', invalid='raise', divide='raise', under='ignore')

def H(q1, q2, p1, p2):
    """
    Função da Hamiltoniana de Hénon-Heiles.
    """
    return 0.5 * (p1**2 + p2**2 + q1**2 + q2**2) + q1**2 * q2 - (1/3) * q2**3

5.2. Implementação do método de Euler simplético e algoritmo de Henon¶

Método de Euler¶

O método de Euler Simplético para esta Hamiltoniana separável é um mapa que atualiza as posições e momentos em duas etapas:

  1. Atualização das Posições: $$q_1^{(n+1)}=q_1^{(n)}+\Delta t p_1^{(n)}$$ $$q_2^{(n+1)}=q_2^{(n)}+\Delta t p_2^{(n)}$$

  2. Atualização dos Momentos (usando as novas posições): $$p_1^{(n+1)}=p_1^{(n)}-\Delta t (q_1^{(n+1)} + 2q_1^{(n+1)}q_2^{(n+1)})$$ $$p_2^{(n+1)}=p_2^{(n)}-\Delta t (q_2^{(n+1)} + (q_1^{(n+1)})^2 - (q_2^{(n+1)})^2)$$

In [7]:
def euler_symplectic_step(q1, q2, p1, p2, dt):
    """
    Realiza um passo de integração usando o método de Euler Simplético.
    """
    # Atualiza as posições
    q1_new = q1 + dt * p1
    q2_new = q2 + dt * p2
    
    # Atualiza os momentos usando as novas posições
    p1_new = p1 - dt * (q1_new + 2 * q1_new * q2_new)
    p2_new = p2 - dt * (q2_new + q1_new**2 - q2_new**2)
    
    return q1_new, q2_new, p1_new, p2_new

Algoritmo de Henon¶

O algoritmo de Henon encontra o ponto exato onde a trajetória cruza a seção de Poincaré $q_1 = 0$. O algoritmo de Hénon envolve um segundo passo de integração para encontrar o ponto exato de interseção com a superfície de seção. Para isso, reformulamos as equações de movimento, usando $q_1$ como a variável independente:

$$\frac{dq_2}{dq_1} = \frac{\dot{q_2}}{\dot{q_1}} = \frac{p_2}{p_1}$$ $$\frac{dp_2}{dq_1} = \frac{\dot{p_2}}{\dot{q_1}} = \frac{-(q_2 + q_1^2 - q_2^2)}{p_1}$$

  1. Verifica o cruzamento
    Checa se $q_1$ mudou de sinal (foi de negativo para positivo) e se $p_1 \ge 0$. Isso garante que a trajetória realmente passou pela seção na direção certa.

  2. Calcula o quanto falta em $q_1$
    A distância até a seção é simplesmente
    $$ \Delta q_1 = -q_1^{(n)} . $$

  3. Avança as outras variáveis em função de $q_1$ Em vez de usar o tempo $t$, o algoritmo usa $q_1$ como “passo”.

  4. Retorna o ponto de interseção

  5. O par $(q_2^{*}, p_2^{*})$ é o ponto preciso na seção $q_1 = 0$.

In [8]:
def henon_alg(q1_old, q2_old, p1_old, p2_old, q1, p1):
    # Condição de cruzamento
    if q1_old < 0 and q1 >= 0 and p1 >= 0:
        # Δq1 necessário para alcançar q1=0
        d_q1 = -q1_old

        # Euler em função de q1
        q2_intersec = q2_old + d_q1 * (p2_old / p1_old)
        p2_intersec = p2_old + d_q1 * (-(q2_old + q1_old**2 - q2_old**2) / p1_old)

        return q2_intersec, p2_intersec

5.3. Seção de Poincaré¶

Esta célula contém a lógica principal para a simulação. A função get_poincare_points referencia o integrador de Euler Simplético e o algoritmo de Hénon para detectar e registrar os pontos da seção de Poincaré.

A seção de Poincaré será construída no plano $(q_2, p_2)$ para a superfície de seção $q_1=0$, com a condição adicional de $p_1 \ge 0$.

In [9]:
def get_poincare_points(initial_state, E_target, dt, num_steps):
    # (código da sua função)
    q1, q2, p1, p2 = initial_state

    poincare_points_q2 = []
    poincare_points_p2 = []
    
    # Loop de integração
    i=0
    for n in range(num_steps):
        # Armazena o estado anterior para detecção de cruzeiro
        q1_old, q2_old, p1_old, p2_old = q1, q2, p1, p2
        
        # Passo de integração principal
        try:
            q1, q2, p1, p2 = euler_symplectic_step(q1_old, q2_old, p1_old, p2_old, dt)
        except FloatingPointError:
            i += 1
            break # Aborta essa trajetória 
        
        # Verifica a condição de cruzamento
        hit = henon_alg(q1_old, q2_old, p1_old, p2_old, q1, p1)

        if hit is not None:
            q2_intersec, p2_intersec = hit 
            poincare_points_q2.append(q2_intersec)
            poincare_points_p2.append(p2_intersec)
        
    return poincare_points_q2, poincare_points_p2, i

5.3. Seção de Poincaré para E=0.08333¶

Para construir a seção de Poincaré, precisamos de uma condição inicial que satisfaça a energia total $E=0.08333$. Vamos escolher uma condição inicial com $q_1=0$ para começar, e então usar a equação da Hamiltoniana para encontrar o valor de $p_2$.

$H = E \implies \frac{1}{2}(p_{1}^{2}+p_{2}^{2}+q_{1}^{2}+q_{2}^{2})+q_{1}^{2}q_{2}-\frac{1}{3}q_{2}^{3} = E$

Fixando $q_1=0$, $p_1=0.2$ e $q_2=0$, podemos resolver para $p_2$.

$E = 0.08333$

$0.08333 = \frac{1}{2}(0.2^2 + p_2^2 + 0^2 + 0^2) + 0 - 0$ $0.08333 = \frac{1}{2}(0.04 + p_2^2)$ $0.16666 = 0.04 + p_2^2$ $p_2^2 = 0.12666 \implies p_2 \approx \pm 0.35589$

Vamos usar $p_2 = 0.35589$.

Figura 4 da Ref. [4]

In [ ]:
# Parâmetros da simulação
E_target = 0.08333
dt = 0.01
num_steps = 10**6

# Geração de múltiplas condições iniciais em uma grade
# Aumente o número de pontos para q2_vals e p2_vals para ter mais trajetórias
# Por exemplo, uma grade 30x30 resulta em 900 combinações
q2_vals = np.linspace(-0.5, 0.5, 30)
p2_vals = np.linspace(-0.5, 0.5, 30)
total_trajectories = len(q2_vals) * len(p2_vals)

all_q2_points = []
all_p2_points = []
num_valid_initials = 0
trajectory_counter = 0

print(f"Iniciando simulações para E = {E_target}...")

# O loop principal para gerar as múltiplas trajetórias
for q2_0 in q2_vals:
    for p2_0 in p2_vals:
        trajectory_counter += 1

        # A condição inicial de p1 é calculada a partir da energia total
        # H = 1/2(p1^2 + p2^2 + q1^2 + q2^2) + q1^2*q2 - 1/3*q2^3
        # Para q1=0, a equação da energia se torna E = 1/2(p1^2 + p2^2 + q2^2) - 1/3*q2^3
        # E pode ser reescrita para isolar p1: p1^2 = 2E - (p2^2 + q2^2) + (2/3)q2^3
        p1_sq = 2 * E_target - (p2_0**2 + q2_0**2) + (2/3) * q2_0**3

        if p1_sq >= 0:
            p1_calc = np.sqrt(p1_sq)
            initial_state = (0.0, q2_0, p1_calc, p2_0)

            # print(f"  > E = {E_target}: Trajetória {trajectory_counter}/{total_trajectories}...")
            q2_points, p2_points = get_poincare_points(initial_state, E_target, dt, num_steps)

            all_q2_points.extend(q2_points)
            all_p2_points.extend(p2_points)

            num_valid_initials += 1

# Bloco para plotar o gráfico
if num_valid_initials > 0:
    plt.figure(figsize=(10, 10))
    plt.scatter(all_q2_points, all_p2_points, s=0.01, alpha=0.5)
    
    # Define os limites de visualização fixos
    q2_min, q2_max = -0.6, 0.6
    p2_min, p2_max = -0.6, 0.6
    plt.xlim(q2_min, q2_max)
    plt.ylim(p2_min, p2_max)

    plt.xlabel('$q_2$', fontsize=12)
    plt.ylabel('$p_2$', fontsize=12)
    plt.title(f'Seção de Poincaré (E = {E_target})', fontsize=14)
    plt.grid(True)

    plt.show()

    print(f"Concluído para E = {E_target}! Total de pontos: {len(all_q2_points)}\n")

Iniciando simulações para E = 0.08333...
No description has been provided for this image 
Concluído para E = 0.08333! Total de pontos: 718199

5.4. Seção de Poincaré para E=0.125¶

Figura 5 da Ref. [4]

In [ ]:
# Parâmetros da simulação
E_target = 0.125
dt = 0.01
num_steps = 10**6

# Geração de múltiplas condições iniciais em uma grade
# Aumente o número de pontos para q2_vals e p2_vals para ter mais trajetórias
# Por exemplo, uma grade 30x30 resulta em 900 combinações
q2_vals = np.linspace(-0.5, 0.5, 30)
p2_vals = np.linspace(-0.5, 0.5, 30)
total_trajectories = len(q2_vals) * len(p2_vals)

all_q2_points = []
all_p2_points = []
num_valid_initials = 0
trajectory_counter = 0

print(f"Iniciando simulações para E = {E_target}...")

# O loop principal para gerar as múltiplas trajetórias
for q2_0 in q2_vals:
    for p2_0 in p2_vals:
        trajectory_counter += 1

        # A condição inicial de p1 é calculada a partir da energia total
        # H = 1/2(p1^2 + p2^2 + q1^2 + q2^2) + q1^2*q2 - 1/3*q2^3
        # Para q1=0, a equação da energia se torna E = 1/2(p1^2 + p2^2 + q2^2) - 1/3*q2^3
        # E pode ser reescrita para isolar p1: p1^2 = 2E - (p2^2 + q2^2) + (2/3)q2^3
        p1_sq = 2 * E_target - (p2_0**2 + q2_0**2) + (2/3) * q2_0**3

        if p1_sq >= 0:
            p1_calc = np.sqrt(p1_sq)
            initial_state = (0.0, q2_0, p1_calc, p2_0)

            # print(f"  > E = {E_target}: Trajetória {trajectory_counter}/{total_trajectories}...")
            q2_points, p2_points = get_poincare_points(initial_state, E_target, dt, num_steps)

            all_q2_points.extend(q2_points)
            all_p2_points.extend(p2_points)

            num_valid_initials += 1

# Bloco para plotar o gráfico
if num_valid_initials > 0:
    plt.figure(figsize=(10, 10))
    plt.scatter(all_q2_points, all_p2_points, s=0.01, alpha=0.5)
    
    # Define os limites de visualização fixos
    q2_min, q2_max = -0.6, 0.7
    p2_min, p2_max = -0.6, 0.7
    plt.xlim(q2_min, q2_max)
    plt.ylim(p2_min, p2_max)

    plt.xlabel('$q_2$', fontsize=12)
    plt.ylabel('$p_2$', fontsize=12)
    plt.title(f'Seção de Poincaré (E = {E_target})', fontsize=14)
    plt.grid(True)

    plt.show()

    print(f"Concluído para E = {E_target}! Total de pontos: {len(all_q2_points)}\n")

Iniciando simulações para E = 0.125...
No description has been provided for this image 
Concluído para E = 0.125! Total de pontos: 1009720

5.5. Séries temporais¶

Séries temporais de $q_1$, $p_1$, $q_2$, $q_2$, para E=0.125 de trajetórias periódicas.

In [19]:
import numpy as np
import matplotlib.pyplot as plt

# Hamiltoniana Hénon–Heiles (para monitorar energia, opcional)
def H(q1, q2, p1, p2):
    return 0.5*(p1*p1 + p2*p2 + q1*q1 + q2*q2) + q1*q1*q2 - (1.0/3.0)*q2**3

# Euler simplético (Euler–Cromer): atualiza q com p^n e p com grad V em q^{n+1}
def euler_symplectic_step(q1, q2, p1, p2, dt):
    # q^{n+1}
    q1_new = q1 + dt*p1
    q2_new = q2 + dt*p2
    # grad V(q1,q2) = ( q1 + 2*q1*q2 , q2 + q1^2 - q2^2 )
    p1_new = p1 - dt*(q1_new + 2.0*q1_new*q2_new)
    p2_new = p2 - dt*(q2_new + q1_new**2 - q2_new**2)
    return q1_new, q2_new, p1_new, p2_new

# gera séries temporais
def simulate_time_series(initial_state, dt, num_steps):
    q1, q2, p1, p2 = initial_state
    t = np.arange(num_steps+1)*dt

    q1_series = np.empty(num_steps+1)
    q2_series = np.empty(num_steps+1)
    p1_series = np.empty(num_steps+1)
    p2_series = np.empty(num_steps+1)

    q1_series[0], q2_series[0], p1_series[0], p2_series[0] = q1, q2, p1, p2

    for n in range(num_steps):
        q1, q2, p1, p2 = euler_symplectic_step(q1, q2, p1, p2, dt)
        q1_series[n+1], q2_series[n+1], p1_series[n+1], p2_series[n+1] = q1, q2, p1, p2

    return t, q1_series, p1_series, q2_series, p2_series

def plot_time_series(t, q1, p1, q2, p2, E_label=None):
    fig, axs = plt.subplots(2, 2, figsize=(11, 6), sharex=True)
    axs = axs.ravel()
    axs[0].plot(t, q1, lw=1); axs[0].set_ylabel(r"$q_1$")
    axs[1].plot(t, p1, lw=1); axs[1].set_ylabel(r"$p_1$")
    axs[2].plot(t, q2, lw=1); axs[2].set_ylabel(r"$q_2$")
    axs[3].plot(t, p2, lw=1); axs[3].set_ylabel(r"$p_2$")
    axs[3].set_xlabel("t")
    axs[2].set_xlabel("t")
    title = "Séries temporais (Euler simplético)"
    if E_label is not None: title += f" — E = {E_label}"
    fig.suptitle(title, y=0.98, fontsize=14)
    for ax in axs: ax.grid(True, alpha=0.3)
    plt.tight_layout(); plt.show()
In [20]:
E_1 = 0.125
dt = 1e-3        # menor para séries suaves
num_steps = 200000  # ajuste para cobrir tempo suficiente

q1_0, q2_0 = 0.0, 0.0
p2_0 = np.sqrt(0.21)      # ~0.45826
# p1 a partir da energia (para q1=0)
p1_0 = np.sqrt(2*E_1 - (p2_0**2 + q2_0**2) + (2/3.0)*q2_0**3)

initial_state = (q1_0, q2_0, p1_0, p2_0)

t, q1_t, p1_t, q2_t, p2_t = simulate_time_series(initial_state, dt, num_steps)
plot_time_series(t, q1_t, p1_t, q2_t, p2_t, E_label=E_1)

print("Energia inicial:", H(*initial_state))
print("Energia final  :", H(q1_t[-1], q2_t[-1], p1_t[-1], p2_t[-1]))
No description has been provided for this image 
Energia inicial: 0.125
Energia final  : 0.12504638709312946

5.6. Séries temporais¶

Séries temporais para E=0.125, trajetória caótica.

In [23]:
import numpy as np
import matplotlib.pyplot as plt

def H(q1, q2, p1, p2):
    return 0.5*(p1*p1 + p2*p2 + q1*q1 + q2*q2) + q1*q1*q2 - (1.0/3.0)*q2**3

def euler_symplectic_step(q1, q2, p1, p2, dt):
    q1_new = q1 + dt*p1
    q2_new = q2 + dt*p2
    # grad V(q) = ( q1 + 2*q1*q2 , q2 + q1**2 - q2**2 )
    p1_new = p1 - dt*(q1_new + 2.0*q1_new*q2_new)
    p2_new = p2 - dt*(q2_new + q1_new**2 - q2_new**2)
    return q1_new, q2_new, p1_new, p2_new

def simulate_time_series(initial_state, dt, num_steps):
    q1, q2, p1, p2 = initial_state
    t  = np.arange(num_steps+1)*dt
    q1s = np.empty(num_steps+1); q2s = np.empty(num_steps+1)
    p1s = np.empty(num_steps+1); p2s = np.empty(num_steps+1)
    q1s[0], q2s[0], p1s[0], p2s[0] = q1, q2, p1, p2
    for n in range(num_steps):
        q1, q2, p1, p2 = euler_symplectic_step(q1, q2, p1, p2, dt)
        q1s[n+1], q2s[n+1], p1s[n+1], p2s[n+1] = q1, q2, p1, p2
    return t, q1s, p1s, q2s, p2s
In [24]:
E = 0.16667
q1_0, q2_0 = 0.0, 0.0
p1_0       = 0.2
p2_0       = np.sqrt(0.29334)  # ≈ 0.5416  (compatível com E=0.16667)

dt = 1e-3
num_steps = 300_000

t, q1_t, p1_t, q2_t, p2_t = simulate_time_series((q1_0, q2_0, p1_0, p2_0), dt, num_steps)

fig, axs = plt.subplots(2,2, figsize=(11,6), sharex=True)
axs = axs.ravel()
axs[0].plot(t, q1_t, lw=0.8); axs[0].set_ylabel(r"$q_1$")
axs[1].plot(t, p1_t, lw=0.8); axs[1].set_ylabel(r"$p_1$")
axs[2].plot(t, q2_t, lw=0.8); axs[2].set_ylabel(r"$q_2$"); axs[2].set_xlabel("t")
axs[3].plot(t, p2_t, lw=0.8); axs[3].set_ylabel(r"$p_2$"); axs[3].set_xlabel("t")
fig.suptitle(f"Séries temporais — trajetória caótica (E={E})", y=0.98)
for ax in axs: ax.grid(True, alpha=0.3)
plt.tight_layout(); plt.show()
No description has been provided for this image
In [6]:
E_2 = 0.125
dt = 0.01
num_steps = 10**6

# Condição inicial para E = 0.16667
# 0.16667 = 0.5 * (0.2**2 + p_2**2)
# 0.33334 = 0.04 + p_2**2
# p_2**2 = 0.29334 => p_2 = sqrt(0.29334) ~ 0.5416
q1_0, q2_0, p1_0, p2_0 = 0.0, 0.0, 0.2, np.sqrt(0.29334)

initial_state = (q1_0, q2_0, p1_0, p2_0)
q2_points, p2_points = get_poincare_points(initial_state, E_2, dt, num_steps)

# Plotar os resultados
plt.figure(figsize=(8, 8))
plt.scatter(q2_points, p2_points, s=1)
plt.xlabel('$q_2$', fontsize=12)
plt.ylabel('$p_2$', fontsize=12)
plt.title(f'Seção de Poincaré (E = {E_2})', fontsize=14)
plt.grid(True)
plt.show()

print(f"Número de pontos na seção: {len(q2_points)}")
No description has been provided for this image 
Número de pontos na seção: 1414

5.6. Seção de Poincaré para E=0.16667¶

In [16]:
# Parâmetros da simulação
E_target = 0.16667
dt = 0.01
num_steps = 10**6

# Geração de múltiplas condições iniciais em uma grade
# Aumente o número de pontos para q2_vals e p2_vals para ter mais trajetórias
# Por exemplo, uma grade 30x30 resulta em 900 combinações
q2_vals = np.linspace(-0.5, 0.5, 30)
p2_vals = np.linspace(-0.5, 0.5, 30)
total_trajectories = len(q2_vals) * len(p2_vals)

all_q2_points = []
all_p2_points = []
num_valid_initials = 0
trajectory_counter = 0

print(f"Iniciando simulações para E = {E_target}...")

# O loop principal para gerar as múltiplas trajetórias
i=0
for q2_0 in q2_vals:
    for p2_0 in p2_vals:
        trajectory_counter += 1

        # A condição inicial de p1 é calculada a partir da energia total
        # H = 1/2(p1^2 + p2^2 + q1^2 + q2^2) + q1^2*q2 - 1/3*q2^3
        # Para q1=0, a equação da energia se torna E = 1/2(p1^2 + p2^2 + q2^2) - 1/3*q2^3
        # E pode ser reescrita para isolar p1: p1^2 = 2E - (p2^2 + q2^2) + (2/3)q2^3
        p1_sq = 2 * E_target - (p2_0**2 + q2_0**2) + (2/3) * q2_0**3

        if p1_sq >= 0:
            p1_calc = np.sqrt(p1_sq)
            initial_state = (0.0, q2_0, p1_calc, p2_0)

            # print(f"  > E = {E_target}: Trajetória {trajectory_counter}/{total_trajectories}...")
            q2_points, p2_points, j = get_poincare_points(initial_state, E_target, dt, num_steps)

            i += j
            all_q2_points.extend(q2_points)
            all_p2_points.extend(p2_points)

            num_valid_initials += 1

# Bloco para plotar o gráfico
if num_valid_initials > 0:
    plt.figure(figsize=(10, 10))
    plt.scatter(all_q2_points, all_p2_points, s=0.01, alpha=0.5)
    
    # Define os limites de visualização fixos
    q2_min, q2_max = -0.6, 1.0
    p2_min, p2_max = -0.6, 1.0
    plt.xlim(q2_min, q2_max)
    plt.ylim(p2_min, p2_max)

    plt.xlabel('$q_2$', fontsize=12)
    plt.ylabel('$p_2$', fontsize=12)
    plt.title(f'Seção de Poincaré (E = {E_target})', fontsize=14)
    plt.grid(True)

    plt.show()

    print(f"Concluído para E = {E_target}! Total de pontos: {len(all_q2_points)}\n Total de Trajetórias Abortadas: {i}\n")

Iniciando simulações para E = 0.16667...

/var/folders/g1/_sv3n23d2yq5w04jzglv_r580000gn/T/ipykernel_55173/1170339906.py:10: RuntimeWarning: overflow encountered in scalar multiply
  p1_new = p1 - dt * (q1_new + 2 * q1_new * q2_new)
/var/folders/g1/_sv3n23d2yq5w04jzglv_r580000gn/T/ipykernel_55173/1170339906.py:11: RuntimeWarning: overflow encountered in scalar power
  p2_new = p2 - dt * (q2_new + q1_new**2 - q2_new**2)
/var/folders/g1/_sv3n23d2yq5w04jzglv_r580000gn/T/ipykernel_55173/1170339906.py:11: RuntimeWarning: invalid value encountered in scalar subtract
  p2_new = p2 - dt * (q2_new + q1_new**2 - q2_new**2)
/var/folders/g1/_sv3n23d2yq5w04jzglv_r580000gn/T/ipykernel_55173/1170339906.py:10: RuntimeWarning: invalid value encountered in scalar subtract
  p1_new = p1 - dt * (q1_new + 2 * q1_new * q2_new)
No description has been provided for this image 
Concluído para E = 0.16667! Total de pontos: 979611
 Total de Trajetórias Abortadas: 0

Seção de Poincaré para E variando de 0.0050 a 0.1667¶

Vamos fazer a mesma coisa para cada energia na lista: 0.0050, 0.0070, 0.0089, 0.0100, 0.0117, 0.0133, 0.0148, 0.0167, 0.0183, 0.0188, 0.0217, 0.0250, 0.0267, 0.0300, 0.0333, 0.0350, 0.0383, 0.0417, 0.0433, 0.0467, 0.0500, 0.0517, 0.0550, 0.0583, 0.0600, 0.0633, 0.0667, 0.0683, 0.0717, 0.0750, 0.0767, 0.0800, 0.0833, 0.0850, 0.0883, 0.0917, 0.0933, 0.0967, 0.1000, 0.1017, 0.1050, 0.1083, 0.1100, 0.1133, 0.1167, 0.1183, 0.1217, 0.1250, 0.1267, 0.1300, 0.1333, 0.1350, 0.1383, 0.1417, 0.1433, 0.1467, 0.1500, 0.1517, 0.1550, 0.1583, 0.1600, 0.1633, 0.1667 , variando também as trajetórias em cada energia para 50 valores de q e 50 valores de p gerando 2500 condições iniciais distintas.

In [22]:
E_list = [
    0.0050, 0.0070, 0.0089, 0.0100, 0.0117, 0.0133, 0.0148, 0.0167, 0.0183, 
    0.0188, 0.0217, 0.0250, 0.0267, 0.0300, 0.0333, 0.0350, 0.0383, 0.0417, 
    0.0433, 0.0467, 0.0500, 0.0517, 0.0550, 0.0583, 0.0600, 0.0633, 0.0667, 
    0.0683, 0.0717, 0.0750, 0.0767, 0.0800, 0.0833, 0.0850, 0.0883, 0.0917, 
    0.0933, 0.0967, 0.1000, 0.1017, 0.1050, 0.1083, 0.1100, 0.1133, 0.1167, 
    0.1183, 0.1217, 0.1250, 0.1267, 0.1300, 0.1333, 0.1350, 0.1383, 0.1417, 
    0.1433, 0.1467, 0.1500, 0.1517, 0.1550, 0.1583, 0.1600, 0.1633, 0.1667
]
In [28]:
# Parâmetros da simulação
dt = 0.01
num_steps = 10**6

# Geração de múltiplas condições iniciais em uma grade
# q2_vals = np.linspace(-0.5, 0.5, 50)
# p2_vals = np.linspace(-0.5, 0.5, 50)
# total_trajectories = len(q2_vals) * len(p2_vals)

# Geração de múltiplas condições iniciais em uma grade
q2_vals_grid = np.linspace(-0.5, 0.5, 50)
p2_vals_grid = np.linspace(-0.5, 0.5, 50)
q2_p2_combinations = np.array(np.meshgrid(q2_vals_grid, p2_vals_grid)).T.reshape(-1, 2)
total_combinations = len(q2_p2_combinations)

# Define o número inicial e final de trajetórias
min_trajectories = 350
max_trajectories = 1000
In [30]:
# Cria a pasta 'imagens' para guardar os gráficos de cada seção de Poincaré
output_dir = "imagens"
if not os.path.exists(output_dir):
    os.makedirs(output_dir)
    print(f"Pasta '{output_dir}' criada.")

# Define limites fixos para os eixos X e Y para todos os gráficos
# Isso garante que a escala não mude, o que é ideal para um GIF.
q2_min, q2_max = -0.6, 0.6
p2_min, p2_max = -0.6, 0.6

Pasta 'imagens' criada.
In [ ]:
# A lista de combinações selecionadas, que irá aumentar conforme aumenta a energia
selected_combinations = np.empty((0, 2))
# O conjunto de índices já presentes na simulação
simulated_indices = set()

# Itera sobre a lista de energias com um índice
for i, E_target in enumerate(E_list):

    all_q2_points = []
    all_p2_points = []
    num_valid_initials = 0

    # Calcula o número total de trajetórias que deveriam ser testadas para essa energia
    num_trajectories_target = int(min_trajectories + (max_trajectories - min_trajectories) * (i / (len(E_list) - 1)))
    
    # Adiciona novas trajetórias se a lista atual for menor que o alvo
    trajectories_to_add = num_trajectories_target - len(selected_combinations)
    
    if trajectories_to_add > 0:
        # Pega os índices que ainda não foram usados
        available_indices = np.array(list(set(range(total_combinations)) - simulated_indices))
        
        if len(available_indices) > 0:
            # Seleciona aleatoriamente novos índices
            new_indices = np.random.choice(available_indices, min(trajectories_to_add, len(available_indices)), replace=False)
            
            # Adiciona as novas combinações à lista
            selected_combinations = np.vstack([selected_combinations, q2_p2_combinations[new_indices]])
            
            # Adiciona os novos índices ao conjunto de usados
            simulated_indices.update(new_indices)

    print(f"Iniciando simulações para E = {E_target} com {len(selected_combinations)} trajetórias...")

    trajectory_counter = 0
    for q2_0, p2_0 in selected_combinations:
        trajectory_counter += 1

        p1_sq = 2 * E_target - (p2_0**2 + q2_0**2) + (2/3) * q2_0**3

        if p1_sq >= 0:
            p1_calc = np.sqrt(p1_sq)
            initial_state = (0.0, q2_0, p1_calc, p2_0)

            # print(f"  > E = {E_target}: Trajetória {trajectory_counter}/{len(selected_combinations)}...")
            q2_points, p2_points = get_poincare_points(initial_state, E_target, dt, num_steps)

            all_q2_points.extend(q2_points)
            all_p2_points.extend(p2_points)

            num_valid_initials += 1

    if num_valid_initials > 0:
        plt.figure(figsize=(10, 10))
        plt.scatter(all_q2_points, all_p2_points, s=0.02, alpha=0.5)
        
        # Define os limites de visualização fixos
        plt.xlim(q2_min, q2_max)
        plt.ylim(p2_min, p2_max)

        plt.xlabel('$q_2$', fontsize=12)
        plt.ylabel('$p_2$', fontsize=12)
        plt.title(f'Seção de Poincaré (E = {E_target})', fontsize=14)
        plt.grid(True)
        
        filename = f'simulacao_{E_target:.4f}'.replace('.', '_')
        filepath = os.path.join(output_dir, filename)
        plt.savefig(filepath, bbox_inches='tight')
        
        print(f"Gráfico salvo como '{filepath}'")
        plt.close() # Fecha a figura para não sobrecarregar a memória
        # print(f"Concluído para E = {E_target}! Total de pontos: {len(all_q2_points)}\n")

    else:
        print(f"Nenhuma trajetória válida encontrada para E = {E_target}. Aumente o intervalo da grade ou a energia.\n")

Gif para as simulações¶

In [3]:
import glob

from PIL import Image as PILImage
from IPython.display import Image as DispImage, display

# 1) Gerar o GIF (se já tiver o .gif, pule esta parte)
frames = [Image.open(p) for p in sorted(glob.glob("imagens/*.png"))]
frames[0].save(
    "animacao.gif",
    save_all=True,
    append_images=frames[1:],
    duration=250,   # ms por fram
    loop=0,         # 0 = loop infinito
    disposal=2      # ajuda a evitar “fantasmas” de frames
)

# 2) Exibir EMBUTINDO no output (vai em base64 no HTML exportado)
# display(IPImage(filename="animacao.gif", format="gif", embed=True))
display(DispImage(filename="animacao.gif", format="gif", embed=True))

<IPython.core.display.Image object>
In [5]:
import base64
from IPython.display import HTML
from pathlib import Path

# Lê e converte para base64
data = Path("animacao.gif").read_bytes()
b64 = base64.b64encode(data).decode("ascii")

# Embute direto no output (vai sair no HTML final)
HTML(f'<img src="data:image/gif;base64,{b64}">')
Out[5]:
No description has been provided for this image
In [ ]: